Simulation PE dependent LR

Pre set-up

Import libraries

library(reshape2)
library(ggplot2)
library(plotly)
library(plyr)
library(Rfast)
library(data.table)
library(knitr)
library(rstudioapi)

Import own functions

# Get path current script is in
source_path = dirname(rstudioapi::getSourceEditorContext()$path)
# Source functions required for this script
source(file.path(source_path, 'Beta_pseudo_sim.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Gaussian_pseudo_sim.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Create_design.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Create_miniblock.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'PE_dep_LR_choice_model.R', fsep = .Platform$file.sep))

Parameter set-up

Set up reward space

# Reward lower and upper boundary
reward_space_lb = 1
reward_space_ub = 100
# Set up space for plotting
reward_space_length = 1000
reward_space = seq(reward_space_lb, 
                   reward_space_ub,
                   length=reward_space_length)

Gaussian parameters and distribution

# Parameters
gaussian_mean = (100 * 1/6) * 3
gaussian_sd = (100 * 1/6) / 3
# Distribution
gaussian_densities = dnorm(reward_space, gaussian_mean, gaussian_sd)
# Normalize density
gaussian_densities = scale(gaussian_densities) - min(scale(gaussian_densities))

Beta parameters and distributions

# Set up "skewdness" parameter for easier change of slope
beta_skewedness = 5
# Parameters for lower end beta and upper end beta
beta_a_le = beta_skewedness/3
beta_b_le = 2*beta_skewedness/3
beta_a_ue = 2*beta_skewedness/3
beta_b_ue = beta_skewedness/3
# Create distributions
beta_densities_le = dbeta(seq(0,1,length=reward_space_length), beta_a_le,beta_b_le)
beta_densities_ue = dbeta(seq(0,1,length=reward_space_length), beta_a_ue,beta_b_ue)
# Normalize densities
beta_densities_le = scale(beta_densities_le) - min(scale(beta_densities_le))
beta_densities_ue = scale(beta_densities_ue) - min(scale(beta_densities_le))

Set rare event threshold

# Set rare events to be 20% of the sampled data
rare_lim = 0.2

---

Control plots

Plot distributions

# Form data frame from density values
data_densities = data.frame(reward_space, 
                            gaussian_densities,
                            beta_densities_le,
                            beta_densities_ue)
# Bring data frame into long format
data_densities = melt(data_densities, id.vars='reward_space')
# Disribution plot
p_dist = ggplot(data=data_densities, aes(x=reward_space, y=value, group=as.factor(variable),
                                    color=as.factor(variable))) +
  geom_line(size=1) +
  scale_color_brewer(palette='Accent') +
  scale_x_continuous(limits = c(1,100)) +
  #stat_summary(fun.data = 'mean', geom = 'vline') +
  geom_vline(xintercept=c(1/3, 1/2, 2/3)*100, linetype='dashed', alpha=0.5, size=0.1) +
  labs(title='Density functions of used distributions',
       x='Outcome',
       y='Normalized density',
       color='Distributions') +
  theme(plot.title=element_text(size=15, face='bold', hjust=0.5),
        axis.title=element_text(size=12))
p_dist

Test sampling with histograms

# Define parameters
n_samples = 500
data_sampling = data.frame(c(1:n_samples))
colnames(data_sampling) = 'Sample'
# Fill data frame with samples
data_sampling$stim_1 = Beta_pseudo_sim(n_samples, beta_a_le, beta_b_le, 'b_le',
                                       reward_space_lb, reward_space_ub)$outcome
data_sampling$stim_2 = Gaussian_pseudo_sim(n_samples, gaussian_mean, gaussian_sd, 'g_an',
                                           reward_space_lb, reward_space_ub)$outcome
data_sampling$stim_3 = Beta_pseudo_sim(n_samples, beta_a_ue, beta_b_ue, 'b_ue',
                                       reward_space_lb, reward_space_ub)$outcome
data_sampling = melt(data_sampling, id.vars='Sample')
# Plot histograms
p_sampling = ggplot(data=data_sampling, aes(x=value, fill=variable)) +
  scale_fill_brewer(palette = 'Accent') +
  geom_histogram(binwidth = 1) +
  facet_wrap(~variable)
ggplotly(p_sampling)

---

Simulating data for the choice model